hvCpG_BayesAlgo_exploration

knitr::opts_chunk$set(echo = TRUE)
invisible(lapply(
  c("dplyr", "data.table", "matrixStats", "ggplot2", "reshape2",
    "ggrepel", "pROC", "stringr", "tibble", "viridis", "ggridges"),
  function(pkg) suppressPackageStartupMessages(library(pkg, character.only = TRUE))
))

Data prep

##########################################################################################
## LSHTM: Maria’s data that she used for the hvCpG paper can be accessed from ing-p5 here:
# /mnt/old_user_accounts/p3/maria/PhD/Data/datasets/
#   
# Her hvCpG scripts are here:
# /mnt/old_user_accounts/p3/maria/PhD/Projects/hvCpGs/Scripts/

## Load results of Maria for hvCpG
Marias_hvCpGs <- read.csv("~/2024_hvCpG/03_prepDatasetsMaria/Derakhshan2022_ST5_hvCpG.txt")

## Recreate Maria's datasets
folder1 <- normalizePath("/mnt/old_user_accounts/p3/maria/PhD/Data/datasets/GEO/BMIQ + 10 PCs + age + sex OUTLIERS REMOVED/")
rds_files1 <- list.files(path = folder1, pattern = "\\.RDS$", full.names = TRUE)

# Read all .rds files into a list and convert each to a matrix
rds_list_mat1 <- lapply(rds_files1, function(file) {as.matrix(readRDS(file))})
## Name the list elements by file names (without extensions)
names(rds_list_mat1) <- gsub("\\.RDS$", "", basename(rds_files1))

rds_list_mat2 <- readRDS("/mnt/old_user_accounts/p3/maria/PhD/Data/datasets/TCGA/TCGA_BMIQ_age_sex_PC_adjusted_OUTLIERS_REMOVED_round2.RDS")
rds_list_mat2 <- lapply(rds_list_mat2, as.matrix)

rds_list_mat <- c(rds_list_mat1, rds_list_mat2) ; rm(rds_list_mat1, rds_list_mat2)

# Keep only the array background (covered in 15 datasets out of 30)
all_cpgs <- unlist(lapply(rds_list_mat, rownames))
cpg_counts <- table(all_cpgs)
common_cpgs <- names(cpg_counts[cpg_counts >= 15])

rds_list_mat <- lapply(rds_list_mat, function(mat) {
  mat[rownames(mat) %in% common_cpgs, ]
})

## Load mQTL-matched controls
cistrans_GoDMC_hvCpG_matched_control <- 
  read.table("~/2024_hvCpG/03_prepDatasetsMaria/cistrans_GoDMC_hvCpG_matched_control.txt", header = T)

## Load cpgnames (all dpg covered after filtration)
cpgnames <- read.table("~/arraysh5files/sorted_common_cpgs.txt")$V1

## Check if I'm consistent
table(cpgnames %in% common_cpgs) ## PERFECT

## 
##   TRUE 
## 406036

sub_cistrans_GoDMC_hvCpG_matched_control <- cistrans_GoDMC_hvCpG_matched_control[
  cistrans_GoDMC_hvCpG_matched_control$hvCpG_name %in% cpgnames &
    cistrans_GoDMC_hvCpG_matched_control$controlCpG_name %in% cpgnames,]

rm(common_cpgs, cpg_counts,all_cpgs, folder1, rds_files1, cistrans_GoDMC_hvCpG_matched_control)

Reproduce Maria’s results

Maria’s paper: We defined an hvCpG in the following way:

1. in 65% of datasets in which the CpG is covered (following quality control), it has methylation variance in the top 5% of all (non-removed) CpGs.
2. is covered in at least 15 of the 30 datasets.

## Within each dataset, calculate the CpGs variance, and keep the top 5%
top5pcvar <- lapply(rds_list_mat, function(mat) {
  # Step 1: Calculate row variances
  row_variances = rowVars(mat, na.rm = T)
  df = data.frame(var=row_variances, cpg=names(row_variances))
  # Step 2
  top = top_n(df, as.integer(0.05*nrow(df)), var) ## slightly different than quantile cause keep ties
  return(top)
})

## CpGs in the top 5% variance:
cpg_counts_top <- table(unlist(lapply(top5pcvar, rownames))) %>%
  data.frame() %>% dplyr::rename("cpgs"="Var1", "all_cpgs_top5pc"="Freq")

## all covered CpGs:
cpg_counts <- table(unlist(lapply(rds_list_mat, rownames))) %>%
  data.frame() %>% dplyr::rename("cpgs"="Var1", "all_cpgs"="Freq")

## Both:
cpg_counts_full <- merge(cpg_counts, cpg_counts_top, all.x = T)
rm(cpg_counts, cpg_counts_top)

## in 65% of datasets in which the CpG is covered (following quality control), it has methylation variance in the top 5% of all (non-removed) CpGs:
## rounding, to mimick Maria's approach
hvCpGs_repro_maria <- na.omit(cpg_counts_full[round(cpg_counts_full$all_cpgs_top5pc/cpg_counts_full$all_cpgs, 2) >= 0.65,])

Maria detected 4143 hvCpGs. I detected 4377 hvCpGs. We both have 4143 hvCpGs in common.

Identify hvCpGs based on a sd multiplicative factor lambda

The proportion of CpGs than Maria found hvCpGs is p(hvCpG) = 1.02%.

Calculate lambda

Within each dataset k, calculate the median sd of all CpG j

all_sd_jk <- sapply(rds_list_mat, function(k){
  
  get_sd_k <- function(k){
    ## Calculate a vector of the row (=per CpG j) sd
    return(rowSds(k, na.rm = T))
  }
  ## Return a vector of sds, in a list for each dataset 
  return(get_sd_k(k))
})

## Plot:
df_long <- reshape2::melt(all_sd_jk, variable.name = "Vector", value.name = "SDs")

## Plot distributions of all vectors on the same graph
ggplot(df_long, aes(x = SDs, color = L1)) +
  geom_density() +
  labs(title = "Distribution of SDs accross datasets",
       x = "SDs",
       y = "Density") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_x_sqrt()

## Calculate lambda per dataset
lambdas = sapply(all_sd_jk, function(x){quantile(x, 0.95, na.rm=T)/median(x, na.rm=T)})
names(lambdas) <- gsub(".95%", "", names(lambdas))

# Histogram with kernel density
ggplot(data.frame(lambda=lambdas, dataset=names(lambdas)),
       aes(x = lambdas)) +
  geom_histogram(aes(y = after_stat(density)),
                 colour = 1, fill = "white", binwidth = .1) +
  geom_density(lwd = 1, colour = 4,
               fill = 4, alpha = 0.25)+
  theme_minimal() +
  geom_vline(xintercept = median(lambdas), col = "red")

median(lambdas) #1.878303

## [1] 4.192186

Which datasets have a higher lambda, and why?

df = data.frame(nind=sapply(rds_list_mat, ncol),
                dataset=names(rds_list_mat),
                lambda=lambdas, 
                tissue = sapply(strsplit(names(rds_list_mat), "_"), `[`, 1),
                ethnicity=sapply(strsplit(names(rds_list_mat), "_"), `[`, 2)) %>%
  mutate(dataset=ifelse(dataset %in% c("Blood_Cauc", "Blood_Hisp"), "Blood_Cauc_Hisp", dataset)) %>% 
  mutate(dataset=ifelse(dataset %in% c("Blood_Mexican", "Blood_PuertoRican "), "Blood_Mex_PuertoRican ", dataset)) %>% 
  mutate(dataset=ifelse(dataset %in% c("CD4+_Estonian", "CD8+_Estonian"), "CD4+_CD8+_Estonian", dataset)) %>% 
  mutate(dataset=ifelse(dataset %in% c("Saliva_Hisp", "Saliva_Cauc"), "Saliva_Hisp_Cauc", dataset))

ggplot(data = df,
       aes(x=lambda, y=nind))+
  geom_smooth(method = "lm")+
  geom_point()+
  geom_label_repel(aes(label = dataset, fill=dataset), size= 2, alpha=.8, max.overlaps = 25)+
  theme_bw()+theme(legend.position = "none")+
  scale_x_log10()

## `geom_smooth()` using formula = 'y ~ x'

## Emphasize the outliers
df_long$col = "grey"
df_long$col[df_long$L1 %in% "BulkFrontalCortex"] <- "red"

ggplot(df_long, aes(x = SDs, group = L1, fill = col)) +
  geom_density(data = df_long[!df_long$L1 %in% "BulkFrontalCortex",], alpha = .5) +
  geom_density(data = df_long[df_long$L1 %in% "BulkFrontalCortex",], alpha = .6) +
  scale_fill_manual(values = c("grey", "red"))+
  labs(title = "Distribution of SDs accross datasets",subtitle = "Red=BulkFrontalCortex",
       x = "SDs",
       y = "Density") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_x_sqrt()

Maximum likelihood analysis

The equation is:

\[ \log\left(P(M_j)\right) = \sum_{i=1}^{n} \log\left( \sum_{Z_j=0}^{1} \left( \sum_{Z_{j,k}=0}^{1} P(M_{i,j} \mid Z_{j,k}) \times P(Z_{j,k} \mid Z_j) \right) \times P(Z_j) \right) \]

Examine weird values

The transformation lead to weird values:

test <- readRDS("/mnt/old_user_accounts/p3/maria/PhD/Data/datasets/GEO/Raw_cleaned_beta_matrices_GEO/Blood_Hisp")
test[rownames(test) %in% "cg23089912",5:10]

##            GSM1870975 GSM1870977 GSM1870979 GSM1870986 GSM1870989 GSM1870992
## cg23089912  0.8486523  0.8883389  0.8610275          0  0.8572614  0.8153023

test <- readRDS("/mnt/old_user_accounts/p3/maria/PhD/Data/datasets/GEO/BMIQ + 10 PCs + age + sex OUTLIERS REMOVED/Blood_Hisp.RDS")
test[rownames(test) %in% "cg23089912",5:10]

##             GSM1870975  GSM1870977 GSM1870979   GSM1870986  GSM1870989
## cg23089912 0.001379438 0.002837506  0.0086327 1.124546e-17 0.002056392
##              GSM1870992
## cg23089912 0.0001942933

# GSM1870986
# 0.00000000000000001124546 --> give extreme value in scaling, and p dnorm=0

The zero was weirdly transformed. Is that expected?

Optimisation over multiple CpGs

Run on LSHTM server (need to be outside of RStudio for parallelising and h5 files access): source(“S01_run_hvCpGdetection_Marias.R”) ## for different datasets

“Nelder-Mead” or “L-BFGS-B” optimisation method?

## Load data & functions specific to Maria's arrays (to do before)                                                 
system.time(source("../03_prepDatasetsMaria/S02.formatArraysforR.R"))

## ✅ bladder: median_sd=0.0097, lambda=4.3565
## ✅ Blood_Cauc: median_sd=0.0173, lambda=5.0770
## ✅ Blood_Gamb: median_sd=0.0098, lambda=5.1092
## ✅ Blood_Hisp: median_sd=0.0140, lambda=4.8656
## ✅ Blood_Japan: median_sd=0.0151, lambda=3.2459
## ✅ Blood_Mexican: median_sd=0.0132, lambda=3.2132
## ✅ Blood_PuertoRican: median_sd=0.0138, lambda=3.2670
## ✅ breast: median_sd=0.0158, lambda=3.8669
## ✅ Buccals_Cauc: median_sd=0.0177, lambda=3.2204
## ✅ Buccals_Sing_9mo: median_sd=0.0150, lambda=4.1840
## ✅ BulkFrontalCortex: median_sd=0.0174, lambda=4.8999
## ✅ CD4+_Estonian: median_sd=0.0162, lambda=2.6741
## ✅ CD8+_Estonian: median_sd=0.0185, lambda=2.6285
## ✅ colon: median_sd=0.0126, lambda=3.8457
## ✅ CordBlood: median_sd=0.0120, lambda=2.9670
## ✅ head and neck: median_sd=0.0151, lambda=4.3256
## ✅ kidney: median_sd=0.0113, lambda=4.4279
## ✅ liver: median_sd=0.0131, lambda=4.3550
## ✅ lung: median_sd=0.0107, lambda=4.4855
## ✅ PBMC_AfricanAmerican: median_sd=0.0129, lambda=3.4579
## ✅ Placenta: median_sd=0.0223, lambda=4.9945
## ✅ prostate: median_sd=0.0137, lambda=4.9652
## ✅ Saliva_Born: median_sd=0.0117, lambda=4.2006
## ✅ Saliva_Cauc: median_sd=0.0167, lambda=3.3665
## ✅ Saliva_Hisp: median_sd=0.0148, lambda=3.2444
## ✅ Saliva_Ken: median_sd=0.0114, lambda=4.5440
## ✅ Skin_UKTwin: median_sd=0.0167, lambda=2.5167
## ✅ SuperiorTemporalGyrus: median_sd=0.0150, lambda=2.3478
## ✅ thyroid: median_sd=0.0111, lambda=5.0079
## ✅ uterus: median_sd=0.0109, lambda=5.4185
## List of 30
##  $ bladder              : num 0.00973
##  $ Blood_Cauc           : num 0.0173
##  $ Blood_Gamb           : num 0.00985
##  $ Blood_Hisp           : num 0.014
##  $ Blood_Japan          : num 0.0151
##  $ Blood_Mexican        : num 0.0132
##  $ Blood_PuertoRican    : num 0.0138
##  $ breast               : num 0.0158
##  $ Buccals_Cauc         : num 0.0177
##  $ Buccals_Sing_9mo     : num 0.015
##  $ BulkFrontalCortex    : num 0.0174
##  $ CD4+_Estonian        : num 0.0162
##  $ CD8+_Estonian        : num 0.0185
##  $ colon                : num 0.0126
##  $ CordBlood            : num 0.012
##  $ head and neck        : num 0.0151
##  $ kidney               : num 0.0113
##  $ liver                : num 0.0131
##  $ lung                 : num 0.0107
##  $ PBMC_AfricanAmerican : num 0.0129
##  $ Placenta             : num 0.0223
##  $ prostate             : num 0.0137
##  $ Saliva_Born          : num 0.0117
##  $ Saliva_Cauc          : num 0.0167
##  $ Saliva_Hisp          : num 0.0148
##  $ Saliva_Ken           : num 0.0114
##  $ Skin_UKTwin          : num 0.0167
##  $ SuperiorTemporalGyrus: num 0.015
##  $ thyroid              : num 0.0111
##  $ uterus               : num 0.0109
## List of 30
##  $ bladder              : num 4.36
##  $ Blood_Cauc           : num 5.08
##  $ Blood_Gamb           : num 5.11
##  $ Blood_Hisp           : num 4.87
##  $ Blood_Japan          : num 3.25
##  $ Blood_Mexican        : num 3.21
##  $ Blood_PuertoRican    : num 3.27
##  $ breast               : num 3.87
##  $ Buccals_Cauc         : num 3.22
##  $ Buccals_Sing_9mo     : num 4.18
##  $ BulkFrontalCortex    : num 4.9
##  $ CD4+_Estonian        : num 2.67
##  $ CD8+_Estonian        : num 2.63
##  $ colon                : num 3.85
##  $ CordBlood            : num 2.97
##  $ head and neck        : num 4.33
##  $ kidney               : num 4.43
##  $ liver                : num 4.36
##  $ lung                 : num 4.49
##  $ PBMC_AfricanAmerican : num 3.46
##  $ Placenta             : num 4.99
##  $ prostate             : num 4.97
##  $ Saliva_Born          : num 4.2
##  $ Saliva_Cauc          : num 3.37
##  $ Saliva_Hisp          : num 3.24
##  $ Saliva_Ken           : num 4.54
##  $ Skin_UKTwin          : num 2.52
##  $ SuperiorTemporalGyrus: num 2.35
##  $ thyroid              : num 5.01
##  $ uterus               : num 5.42

##    user  system elapsed 
##   4.204   0.060   4.263

## Load algorithm (30sec)                                                                                                                                                                                                                 
system.time(source("../05_hvCpGalgorithm/hvCpG_algorithm_detection_v3.R"))

##    user  system elapsed 
##   0.085   0.004   0.088

## Load full cpg names
cpg_names <- readLines("~/arraysh5files/sorted_common_cpgs.txt")

length(cpg_names); head(cpg_names)

## [1] 406036

## [1] "cg00000029" "cg00000108" "cg00000109" "cg00000165" "cg00000236"
## [6] "cg00000289"

## [1] 406036
# [1] "cg00000029" "cg00000108" "cg00000109" "cg00000165" "cg00000236" "cg00000289"

## Load hvCpG of Maria and matching mQTL
cistrans_GoDMC_hvCpG_matched_control <- read.table("../03_prepDatasetsMaria/cistrans_GoDMC_hvCpG_matched_control.txt", header = T)

## In which positions are they?
cpg_names[cpg_names %in% cistrans_GoDMC_hvCpG_matched_control$hvCpG_name] %>% length

## [1] 3644

cpg_names[cpg_names %in% cistrans_GoDMC_hvCpG_matched_control$controlCpG_name]  %>% length

## [1] 3453

sub_cistrans_GoDMC_hvCpG_matched_control <- cistrans_GoDMC_hvCpG_matched_control[
  cistrans_GoDMC_hvCpG_matched_control$hvCpG_name %in% cpg_names &
    cistrans_GoDMC_hvCpG_matched_control$controlCpG_name %in% cpg_names,]

## Positions of targets in cpg_names:
match(sub_cistrans_GoDMC_hvCpG_matched_control$hvCpG_name, cpg_names)

##    [1] 286867  11867  47291 360391 132127 149298 380388  67889 247489 262008
##   [11] 280695 387017   4687  44190  97721 114305 154515 162750 269223 285461
##   [21] 401478  51312  78828  94696 168965 173956 188127 190187 202999 209910
##   [31] 213907 247886 265860 265989 292116 305761 367254 391684 396703 402077
##   [41] 379436  11622  26433  78472 162858 239725 323051 391687   4233  11413
##   [51]  31005  34621  48784  52814  55043  63806  64760  91836 110035 111472
##   [61] 115392 127405 150991 153593 156804 165464 167602 168844 169111 174663
##   [71] 186965 199384 217451 232599 253892 256746 272225 274437 277201 277801
##   [81] 288463 290707 303017 320062 337615 341652 345045 354380 361269 361649
##   [91] 362418 366459 396602   3987   9483  15511  23464  26431  34428  46953
##  [101]  48747  62797  86585 112252 123529 136471 141256 153053 156899 163733
##  [111] 166007 174727 203031 206930 222348 224561 243967 260485 269287 284375
##  [121] 289404 293012 322092 322899 332752 340562 344085 357125 358220 363955
##  [131] 380589 385824 392169 398302   2258  11392  30844  39649  40490  44848
##  [141]  53546  55201  55729  87163  89354  97204 110421 117356 122433 130635
##  [151] 131809 132163 135077 139778 141961 147938 154664 172583 173001 185333
##  [161] 191979 193205 207127 213500 223005 242347 242825 244255 244403 247702
##  [171] 249079 281610 285652 294478 304392 306092 313876 326306 338071 344691
##  [181] 369793 370544 375852 379637 383293 390851 393843 405996  35477  39137
##  [191]  53460  55483  65044  74218  89281  94251  94496  97123 111668 116837
##  [201] 142568 143560 148417 166805 217753 224934 225316 228497 232276 233090
##  [211] 240721 251069 267102 269069 283829 285789 308806 310134 318549 320026
##  [221] 330572 332603 348649 364547 385775 388580  18026  47826  81976  82873
##  [231] 141294 166625 178353 199380 209131 250835 289161 384622 386081 399464
##  [241]   1231  13389  16157  16898  22548  25948  25992  30833  31218  31699
##  [251]  32502  32688  33321  34583  36699  38507  45243  45967  47905  48544
##  [261]  53479  57297  59068  59883  61386  64215  68452  70446  72733  77855
##  [271]  87052  88990  97010 100472 101047 105594 110579 110632 112389 122770
##  [281] 123111 123199 123740 125674 141037 145537 148574 149261 151292 153444
##  [291] 154404 155610 156806 159544 161832 162254 162694 163749 169655 171464
##  [301] 177528 187485 193326 194993 195581 198287 200675 201846 206990 208678
##  [311] 210442 217026 224059 228575 229409 230010 234745 236747 238519 243665
##  [321] 245978 247472 254710 255156 255178 256673 258135 259874 260219 265384
##  [331] 268731 271630 274477 278738 280487 280617 284977 285740 287332 294215
##  [341] 294728 295893 297554 299130 300623 301990 302953 305395 306307 306751
##  [351] 306886 306888 309967 310553 314259 314843 322168 329145 330883 335267
##  [361] 337259 337331 337510 339086 339998 341094 343042 343476 345129 348135
##  [371] 350257 350801 359957 363786 364740 365129 366690 368046 368760 370075
##  [381] 370516 370908 374272 374822 375727 381849 381950 383595 384447 385161
##  [391] 386557 386577 387781 390499 392222 392425 398777 398842 399803 402852
##  [401] 403106   1051   1319   2534   2642   7988  10707  12517  14598  15134
##  [411]  19991  22610  27589  30335  30709  31938  33822  37543  38843  39179
##  [421]  44037  50529  51344  52421  55591  56745  58849  61372  66972  67096
##  [431]  69182  71005  73576  74939  76316  77430  87303  87549  89265  95002
##  [441]  95305  96819  97114  98710 102028 103088 103643 103966 104789 105112
##  [451] 105647 105716 119554 121915 123222 127339 129498 133248 135914 137290
##  [461] 140000 141344 141656 143295 146776 147606 151132 153162 158087 159428
##  [471] 159788 161748 163740 168832 169485 171866 175798 177895 181429 183015
##  [481] 183751 186765 193318 196358 199890 200578 208209 216520 217103 225711
##  [491] 228738 231479 231883 234283 234475 236737 236901 237611 237766 238190
##  [501] 238976 240300 241231 243022 243045 250772 251574 252062 252413 255470
##  [511] 255894 257362 261537 264212 265978 268941 270020 270111 272552 274594
##  [521] 276112 276485 282312 284724 284978 286002 287754 291403 291615 292076
##  [531] 293910 295220 296463 297448 299030 304387 308666 310968 312935 314718
##  [541] 316855 317664 321568 325350 331124 332868 333699 337061 337765 342902
##  [551] 345533 346020 346481 347753 347889 348803 351360 354720 355181 356210
##  [561] 356489 365287 365661 368808 368840 370813 372952 373292 376232 380386
##  [571] 381147 381583 385644 390894 390948 396749 401627 402395 402774 404528
##  [581] 404669   1133  34999  38259  38309  47246  57616  59458  62395  69731
##  [591]  77410  81391  88587 100629 104190 112733 116587 122715 123316 127802
##  [601] 135386 142285 145575 147455 154356 160904 201303 207903 255886 269072
##  [611] 277221 284922 295805 297161 304232 312728 319658 319700 331180 333257
##  [621] 345012 347287 356032 357465 365418 370844 380232 390745 397241   1152
##  [631]   1642   2037   2963   4671   4727   4875   4952   5152   5669   6008
##  [641]   8180   8474   8643   8886   8929  10918  13225  14027  14136  14450
##  [651]  14594  14871  15846  15962  16136  16555  16598  18046  18592  19316
##  [661]  19460  19527  21282  21500  21879  22006  22702  23214  23723  27860
##  [671]  27992  30022  30795  31461  32153  32360  34009  35624  36418  37064
##  [681]  37653  38191  39055  39400  39822  41133  41461  42096  43364  43430
##  [691]  44062  44167  44561  44647  45465  46023  46657  47301  47567  49236
##  [701]  50868  53520  54107  55592  56264  56702  56788  56895  57395  57654
##  [711]  57820  57995  58461  58740  59041  59285  60071  61530  62222  64174
##  [721]  65319  66733  67152  67699  68263  68874  69688  70372  74365  76698
##  [731]  76854  79428  79686  79775  80873  81742  82334  82728  85459  86869
##  [741]  86892  87725  88647  88834  89266  89449  90023  90405  90549  90602
##  [751]  90835  91236  91363  91377  91973  91977  92277  92601  92957  93045
##  [761]  93851  93959  95295  95334  95375  96564  96974  97182  97746  98980
##  [771]  99518  99818  99916 100223 100475 100824 102525 102584 102646 103880
##  [781] 104589 105243 105328 105547 105834 106267 106501 109644 110208 110619
##  [791] 111079 111441 111475 112552 113533 113812 114489 115638 115805 116344
##  [801] 116619 116741 116753 117075 119149 120016 122354 122869 123148 124710
##  [811] 125278 125355 126053 126286 126638 126651 127370 128420 128906 130960
##  [821] 131294 134352 134994 136014 136020 136461 136736 136934 137355 137407
##  [831] 137834 138240 139562 141631 141836 141969 142612 142709 143233 143966
##  [841] 144130 145276 146234 146637 146711 147749 148330 148755 149328 150929
##  [851] 151216 151299 152078 152121 152235 153320 153403 153696 154933 156284
##  [861] 156939 157132 157975 158271 158446 159345 159663 159888 160093 160475
##  [871] 161587 162207 164792 166117 167459 167620 167655 168347 168807 169340
##  [881] 169421 169775 170074 170682 172043 172369 172779 173104 174345 174521
##  [891] 174610 175571 175867 175915 177293 177427 178344 178523 178790 179684
##  [901] 180592 181744 182345 182667 182855 183289 183376 183377 183568 183819
##  [911] 184281 184988 185190 185575 186184 186293 187157 187172 188452 189143
##  [921] 189196 189255 189755 190721 190989 192499 193515 196515 196966 198014
##  [931] 198966 199327 199413 199421 201084 201671 201937 201973 202079 202614
##  [941] 202671 203019 203354 204934 205823 205848 207672 209176 211330 213834
##  [951] 214102 214132 214374 214666 215885 215958 217877 218126 218910 219122
##  [961] 219568 220840 221199 222983 223416 224350 224720 225106 225265 225737
##  [971] 226386 229733 231517 231654 232033 232850 233280 233742 234066 234518
##  [981] 236266 236428 236886 238254 238284 238972 240557 240607 241222 241329
##  [991] 241399 242556 243670 244064 246425 247462 247790 248038 249548 249600
## [1001] 249942 250299 250919 252307 253271 253614 255055 255336 255485 255764
## [1011] 256515 259597 260324 260633 261174 262954 263231 263867 264810 264927
## [1021] 266272 267513 268328 268723 269600 269695 270107 270499 270804 270815
## [1031] 271887 272190 274651 274978 276017 276286 276743 277190 277395 277840
## [1041] 277847 277853 278939 279146 279322 279855 280319 282367 283756 284126
## [1051] 285841 287007 287782 288589 288715 288839 291125 291505 291558 292082
## [1061] 292283 292427 293092 293759 293971 294901 295581 296713 298015 298768
## [1071] 299946 299993 300119 300960 301445 301449 301571 301741 301758 302630
## [1081] 304913 305779 305793 306983 307064 307728 307730 308725 309389 309639
## [1091] 310020 310577 310830 310869 312710 313407 314501 314812 314921 315481
## [1101] 315862 317434 317698 317888 318538 321145 321259 322089 322137 323074
## [1111] 323093 323678 324657 325183 326431 326583 326686 326935 327037 327312
## [1121] 328374 328620 329392 331680 332294 334867 336628 337130 337355 337416
## [1131] 337608 337642 339205 340347 340508 342490 343097 344601 345054 346573
## [1141] 346691 347809 348789 349186 349441 349816 350272 350480 351184 351187
## [1151] 351782 352332 353739 353967 354309 354625 354817 357015 357782 358526
## [1161] 359741 360039 360317 360546 361022 363817 363908 364101 364936 364941
## [1171] 365068 365929 366810 367401 368012 370137 372384 372497 373015 373085
## [1181] 374851 375170 375427 375553 375679 375737 376876 376949 377038 377053
## [1191] 378418 379025 379767 379810 379963 381090 383363 384514 385145 385373
## [1201] 386187 388509 388966 389005 390094 391021 391198 391467 391852 392092
## [1211] 393036 393134 394110 395051 395798 395926 395960 396078 396999 397984
## [1221] 398514 399920 399963 400000 400581 400994 402445 402715 403685 404062
## [1231] 404080 404242 404322    455    587   1664   2007   2278   2317   3024
## [1241]   3079   3220   3717   4627   5423   7277   7392   7418   7822   8600
## [1251]   9212   9392   9735  10324  10421  10623  10817  10888  11060  13528
## [1261]  13557  13798  14007  14135  14529  14655  15105  15249  15422  15851
## [1271]  15958  16356  16875  17577  17629  18144  18831  18879  19848  19923
## [1281]  20020  20287  20321  20459  20535  21575  21587  22103  22344  22573
## [1291]  22745  22794  23066  23366  23431  23553  23629  23701  23811  24374
## [1301]  24710  24818  24890  24981  25375  25681  25884  26010  26446  26820
## [1311]  27300  27324  28134  28279  29037  29461  30231  30328  30553  30640
## [1321]  30711  30808  31356  31472  31828  32471  32977  34093  34277  34452
## [1331]  34500  34561  35069  35230  35363  36189  36724  37576  37847  38720
## [1341]  38722  38827  38969  39136  39501  40133  40136  40320  40971  41305
## [1351]  41338  41518  41752  41833  42741  42801  43186  43267  43360  43476
## [1361]  43617  43658  44204  44365  44414  44507  44642  44928  45976  46486
## [1371]  47015  47304  47979  48865  49017  49403  49439  49532  49889  50041
## [1381]  50309  50792  50843  51662  52058  52186  53076  53282  53378  53886
## [1391]  53964  54399  54547  54905  54999  55288  56059  56085  57190  58042
## [1401]  58090  58466  58821  59185  59503  59772  60830  61203  61487  61655
## [1411]  62003  62167  62561  62588  62699  62837  63240  63526  63617  63760
## [1421]  63977  64365  65324  65785  65946  66118  66694  66786  66913  67273
## [1431]  67672  67754  67773  68046  68143  68331  68438  68581  68746  68919
## [1441]  69421  69573  70356  71041  72018  72082  72762  73647  73665  73930
## [1451]  74172  74265  75502  75708  77312  78133  78407  78685  78700  78708
## [1461]  78713  79026  79230  79342  79535  80613  80930  81283  81442  81523
## [1471]  81734  81972  83555  83790  83964  84276  84317  84702  84953  85296
## [1481]  86054  86474  86859  86864  86941  87053  87142  87331  87349  87990
## [1491]  88245  88253  88390  88496  88679  89001  89843  90037  90243  90399
## [1501]  91314  91740  91819  92055  92352  92571  93196  93205  93808  94152
## [1511]  94363  94476  94669  95090  95153  95274  95838  95870  97311  97502
## [1521]  97826  97965  98099  98370  99279  99439  99494  99564  99659  99896
## [1531] 100386 100483 100578 101208 101531 101775 101907 101936 102621 102859
## [1541] 102884 103135 103819 104562 105254 105534 106406 106858 106872 106932
## [1551] 107113 107530 107686 107697 107813 107884 107933 108086 108148 108215
## [1561] 108478 108937 109335 109469 109492 109577 109646 109933 110187 110764
## [1571] 111632 111718 111770 112024 112279 112824 112949 113030 113206 113232
## [1581] 113433 114006 114693 115414 115515 115593 116360 116596 116910 117235
## [1591] 117441 118028 118806 118870 119215 119268 119315 119484 119582 119900
## [1601] 120354 120743 121147 121390 122014 122124 123197 123419 123586 123972
## [1611] 124056 124156 124800 124820 124851 125004 125043 125120 125517 125720
## [1621] 125867 125890 126564 126668 126897 127022 127052 127135 127138 127226
## [1631] 129240 129767 129948 130445 130482 130541 130935 131066 131105 131210
## [1641] 131287 131356 131537 131572 131789 132125 132411 132717 132780 132954
## [1651] 133277 133733 134063 134408 134713 134987 135253 135312 135486 135498
## [1661] 136133 136854 136870 137334 137773 138232 138485 138692 139396 139522
## [1671] 140001 140496 141129 141187 141286 141365 141905 142183 142225 142326
## [1681] 143038 143253 143536 143542 143636 143793 144132 144398 144961 145505
## [1691] 145974 146053 146084 146465 146498 146736 147584 148120 148170 148311
## [1701] 148317 148595 149045 150132 150332 150550 150667 150877 151122 151240
## [1711] 151345 151379 151496 151593 151765 151778 151927 152344 152537 153220
## [1721] 153241 153405 153884 154070 154555 154578 154603 154777 154805 155060
## [1731] 155143 155167 155334 156091 156702 156734 156817 157235 157814 157950
## [1741] 158850 158878 158945 159035 159598 159876 159900 160140 160164 160211
## [1751] 160758 161197 161306 161519 161536 161809 162042 162142 162689 163161
## [1761] 163773 164001 164478 164762 164980 165267 167265 167364 167929 168214
## [1771] 168231 168267 168813 169314 169521 169787 170451 170570 170702 171056
## [1781] 171421 172697 172751 173068 173178 173281 173324 173418 174073 174527
## [1791] 174873 175003 175096 175781 175851 175994 176022 176059 176430 176692
## [1801] 176994 177311 177382 178268 178272 178429 178654 179367 179734 179755
## [1811] 179833 180100 180391 180489 181233 181284 181288 181670 182003 182112
## [1821] 182608 182612 183007 183149 183561 184834 185165 185328 185607 185695
## [1831] 185757 185801 185871 186078 186448 187324 188192 188295 188604 189121
## [1841] 189465 189504 189966 190164 190234 190293 191050 191356 191403 191895
## [1851] 191911 191952 193130 194231 194235 195041 195252 195626 195737 195783
## [1861] 196573 196725 196753 196930 197036 197106 198009 198028 198569 199078
## [1871] 199229 199359 200611 200851 201029 201804 201874 201914 202753 203814
## [1881] 203825 204015 204023 204086 204431 205036 205045 205052 205248 205373
## [1891] 206113 206200 206358 206489 206490 206812 207392 207647 208601 208847
## [1901] 208942 209297 209637 209649 210819 211285 211733 211852 211925 212209
## [1911] 212860 212917 213252 213347 213402 213575 214151 214289 214445 214721
## [1921] 215255 216051 216601 216631 217192 217201 217446 217818 218208 218464
## [1931] 218667 218829 218955 219163 219537 219548 220028 220056 220080 221347
## [1941] 221431 222751 223080 223403 223902 224073 224286 224486 224618 225159
## [1951] 225264 225337 225536 226513 227251 228159 228603 228620 228750 228948
## [1961] 229750 229897 230298 230383 230709 231130 231352 231449 231826 232411
## [1971] 232638 232936 232957 233397 233550 233747 233761 233975 234152 234230
## [1981] 235009 235010 235082 235511 235519 236448 236850 237068 237070 237099
## [1991] 237167 237357 237469 237729 238260 238706 239718 239956 240034 240039
## [2001] 240837 240931 241398 241431 241665 241714 241784 241856 242022 242574
## [2011] 242648 242928 243280 243441 243811 244041 244438 244560 244603 245120
## [2021] 245452 245549 245828 246735 246984 247520 248311 248644 248730 249228
## [2031] 249237 249588 249952 250464 250661 252384 252449 252523 252723 253013
## [2041] 253502 253711 253773 254028 254151 254249 254293 254883 255383 255688
## [2051] 256166 256381 256555 256722 257088 257174 257203 257225 257228 257434
## [2061] 257802 257814 258048 258253 258322 258412 258442 258616 258850 259546
## [2071] 259768 259836 259875 260251 260493 260876 261194 261352 261363 261976
## [2081] 262200 262437 262824 262860 263203 263255 263539 264040 264366 264908
## [2091] 265043 265634 265651 265890 266565 266573 267017 267222 267231 267730
## [2101] 268061 268398 269009 269220 269345 269648 269719 269784 270078 270284
## [2111] 270424 270645 270951 271081 272251 272341 272408 272485 272519 272532
## [2121] 272664 272690 272905 272936 273113 273116 273490 273946 274732 274779
## [2131] 274837 275479 275987 276179 276448 277419 277945 277964 278480 278584
## [2141] 278902 279269 280458 280590 280879 281494 281696 281969 282394 282682
## [2151] 282777 283622 284222 284620 284750 284901 285056 285167 285242 285264
## [2161] 286371 286581 286646 286791 286827 287165 287188 287642 287688 287971
## [2171] 287980 288025 288029 288356 288481 288558 288813 288969 289479 289503
## [2181] 289745 290147 290995 292135 292251 292492 292858 293841 293864 294549
## [2191] 294971 295136 295351 295433 295659 295885 295917 295979 296127 296594
## [2201] 296615 296872 297197 297578 297804 297813 297822 298228 298377 298545
## [2211] 299325 299544 299716 299827 299845 300066 300102 300644 300657 300687
## [2221] 300816 301717 301919 302311 303009 303050 303343 303410 303803 304026
## [2231] 304131 304665 304668 304801 305226 305693 306023 306033 306305 306551
## [2241] 307776 308448 308633 308643 308694 308987 309288 309378 309603 309636
## [2251] 309753 309879 310154 310331 310403 310877 310893 311019 311025 311049
## [2261] 311225 312107 312139 312227 313331 313961 314068 314674 314801 315178
## [2271] 315341 315458 315716 316494 316562 316627 316733 317264 317331 317402
## [2281] 317772 318033 318157 318534 318796 319111 319288 319560 319857 319931
## [2291] 319970 320652 320995 321111 321191 321406 322674 322858 323674 324443
## [2301] 324606 324641 324668 324972 324989 325615 325702 325761 326175 326430
## [2311] 326524 326615 326972 327456 327768 327872 327955 328300 328538 328957
## [2321] 328997 329218 329507 329674 329693 329858 331072 331346 331355 331439
## [2331] 331720 331749 332634 332904 332991 333092 333475 334252 334384 334727
## [2341] 334937 335407 335460 335563 335889 336261 336461 336609 336948 337153
## [2351] 337406 337661 338140 338267 338317 338668 339069 339145 339763 341231
## [2361] 341343 341408 341456 341882 342205 342681 343001 343002 343463 343786
## [2371] 343849 344338 344368 344805 345378 345496 345542 345643 345715 345876
## [2381] 346531 346644 347006 348290 348313 348845 348937 349159 351207 351588
## [2391] 351769 351816 351887 352055 352062 352067 352240 352770 353158 353396
## [2401] 353416 353466 353781 353913 353976 354294 354635 355206 355235 355536
## [2411] 356591 356884 356989 357082 357297 358075 358262 358350 358542 358690
## [2421] 359121 359255 359396 359405 359611 360510 360737 360928 361612 361704
## [2431] 361725 361744 362109 362158 362190 362296 362524 362644 363077 363148
## [2441] 363458 364172 364254 364814 366133 366335 366342 366579 366636 367039
## [2451] 367433 367442 367520 367752 368319 368333 369262 369334 369586 369644
## [2461] 370038 370423 370566 371191 371433 371661 371908 372033 372042 372091
## [2471] 372311 372396 372420 372779 372844 373248 373414 373603 373693 373695
## [2481] 373732 373807 374312 374874 374908 375404 375464 376086 376428 377330
## [2491] 377486 378219 378299 378538 378654 378730 379215 380081 380413 380481
## [2501] 380509 381096 381644 382696 382732 382831 383343 383643 384095 384991
## [2511] 385766 386674 387337 387591 388904 388909 389793 390693 390807 391016
## [2521] 391064 391312 392147 392157 392427 392661 392741 392759 392792 392876
## [2531] 392882 393275 393320 393646 393828 394248 394411 394597 394734 395216
## [2541] 395422 395482 395598 397083 397296 398041 398133 398218 398220 398676
## [2551] 398744 399019 399104 399162 399394 400195 400521 401668 401992 402570
## [2561] 403197 403790 403888 404144 404314 404583 405398 405514 405539     96
## [2571]    452    948   1605   1998   2519   2831   2908   3977   4410   5323
## [2581]   5659   5739   5802   6102   6698   6714   6746   6898   6997   7294
## [2591]   8326   9480   9532   9694  10715  11431  11523  12741  13206  13310
## [2601]  13558  14183  14330  14561  14568  14966  15630  16602  16643  17085
## [2611]  18230  18969  19163  19995  21023  21141  21261  21729  21765  21784
## [2621]  22043  22089  22332  22948  23741  23966  24813  24948  25283  25613
## [2631]  25747  26419  26822  27388  27458  27782  27974  28323  28927  29009
## [2641]  29307  30232  30818  31380  31772  32162  32461  33211  33337  33805
## [2651]  34377  35336  35456  35818  36945  37495  38712  39021  39300  39315
## [2661]  39339  39634  39812  39834  40029  40711  40755  41439  41957  42073
## [2671]  42372  42896  42924  43046  43398  43462  43479  43482  43510  43972
## [2681]  44678  44830  45744  46124  46472  46535  47089  48147  48472  48765
## [2691]  48953  49190  49340  49410  50344  51166  52211  53150  53203  54141
## [2701]  54170  54500  54542  54608  54747  55371  55525  55562  55934  56397
## [2711]  56906  57142  57787  58200  58760  58844  58851  58865  59403  59650
## [2721]  59709  60865  61450  61455  62882  64350  64840  65199  65713  65821
## [2731]  65915  66062  66113  66121  66213  66563  66933  67726  68622  69766
## [2741]  70081  70440  71415  71496  71614  71616  72714  74282  75219  76060
## [2751]  76198  76251  76280  77125  77329  78328  78679  79412  80265  81008
## [2761]  81248  82966  83044  85018  85402  85708  88169  89638  89668  89986
## [2771]  90089  91183  91425  93284  93354  94238  94275  94348  94533  94649
## [2781]  94831  95029  95892  95938  96485  96678  97110  97246  97762  98363
## [2791]  98597  99453 100311 100946 101179 101297 101469 101664 102270 102326
## [2801] 102515 102674 102722 104015 104262 105506 105625 105667 106091 106271
## [2811] 106306 106423 107276 107351 107983 108418 108775 109057 109386 110121
## [2821] 110657 111136 111581 111598 111794 113360 113382 114073 114563 114623
## [2831] 114642 115025 115250 115698 116058 117085 117844 118144 118170 118247
## [2841] 118791 118930 119027 119262 119327 119417 119832 120066 120144 120548
## [2851] 120632 121687 123407 123474 123520 123530 124014 124177 124678 124943
## [2861] 125125 125261 125888 125944 127775 127962 128089 128172 128548 128729
## [2871] 128816 128859 128898 129955 130127 130150 130329 130439 130508 131045
## [2881] 131451 131804 132091 132286 132542 132786 132808 134579 134907 135565
## [2891] 135576 137384 138322 138667 139133 139640 139703 140106 140197 140317
## [2901] 140720 140950 141082 141202 142261 142572 142768 142804 143257 143525
## [2911] 143721 143768 144378 144556 144665 145496 146566 147048 147158 147877
## [2921] 149048 149128 149307 150275 150603 152354 152414 152734 153393 153394
## [2931] 153438 153510 155111 155364 155547 156915 157147 157862 158848 159454
## [2941] 160095 160252 160497 160717 160839 161639 162063 162894 162995 163033
## [2951] 163232 163276 163637 164080 164192 164535 165299 166489 166731 167519
## [2961] 168044 168147 168829 169177 169893 170033 170133 170270 170616 171139
## [2971] 171342 171343 171826 172254 173296 173356 173599 173785 174082 174476
## [2981] 175204 175533 175569 176785 178040 178302 178370 178391 178879 179060
## [2991] 179748 179817 179923 179992 180418 181055 181580 181767 183304 184394
## [3001] 184504 186118 186549 186591 186602 187923 187983 188005 188476 189687
## [3011] 189703 190209 190300 192929 193055 193304 193332 193342 194110 194620
## [3021] 195717 196938 197071 197370 198383 200066 201695 202143 202175 202391
## [3031] 202948 203012 203570 203743 203881 204172 204468 205884 206060 206419
## [3041] 206827 207047 207574 207867 208054 209081 210744 210902 211848 212136
## [3051] 212577 213336 213773 213837 214048 214133 214794 215247 215786 216532
## [3061] 216715 217958 218062 218213 218486 219476 219772 219870 220043 220046
## [3071] 220491 220676 221748 221824 222467 224095 225309 225610 225710 225937
## [3081] 226362 226477 226556 226998 227114 227583 227960 228252 228549 228686
## [3091] 228732 229133 229254 230620 230977 231824 233245 233498 233681 233711
## [3101] 233952 234003 234351 235099 235641 235910 236046 236443 236495 236591
## [3111] 236887 237325 237457 237633 237816 238622 239565 240019 240185 240998
## [3121] 241127 241428 242213 243060 243143 243293 244638 245292 245656 245844
## [3131] 245977 248299 248524 248937 249737 249772 249931 251107 252198 252209
## [3141] 252649 254762 254976 255286 255672 256401 256685 256772 257245 257611
## [3151] 259054 259343 261373 262233 262323 262535 262607 263204 263425 263530
## [3161] 263620 263643 263776 264159 265093 265201 265284 266218 266286 267424
## [3171] 268546 268595 269010 269016 269392 269814 270147 270278 271435 272317
## [3181] 272388 272703 273081 273334 273732 273828 273897 274425 274615 274753
## [3191] 275413 275473 276279 276322 276367 276951 276987 277152 278534 279150
## [3201] 279672 279717 280245 280464 281025 281106 281182 281265 281281 281786
## [3211] 281888 282099 282277 282281 284344 285649 285791 286331 286688 287038
## [3221] 287506 288155 288248 288633 288759 288948 289494 289623 290383 290579
## [3231] 290842 293643 294375 295453 295454 295497 295656 295686 295721 295998
## [3241] 296981 298678 300246 300503 300504 301000 301913 302131 302238 302578
## [3251] 303573 304763 305663 306578 306854 306971 307610 310758 311136 311410
## [3261] 311578 312225 312237 312764 313085 314110 314172 314578 315041 315122
## [3271] 315331 316656 316728 316926 317165 317250 318525 319211 319225 320442
## [3281] 321123 321800 322067 322755 323100 323438 323684 323896 325584 326341
## [3291] 326925 327130 327470 328020 328233 328968 329128 329443 329511 329720
## [3301] 330287 330362 330509 330585 331073 331893 331989 332258 332668 332888
## [3311] 333091 334171 334624 335455 335579 335673 335699 335957 336030 336389
## [3321] 336438 337385 338456 339169 339748 339838 340183 340357 340710 341799
## [3331] 343227 344625 345044 345770 346369 347172 347493 347649 348471 348681
## [3341] 349198 350302 350584 351130 351225 353045 353893 354009 354872 354928
## [3351] 354966 355560 356223 357019 357136 357841 359223 359995 360162 361098
## [3361] 361715 362321 362522 362527 362556 362835 363675 364152 364327 365734
## [3371] 365900 366140 366225 366289 367662 367737 368156 369233 369532 369663
## [3381] 370509 371311 371842 371986 372557 372863 373822 374753 375086 375352
## [3391] 375419 376129 376181 376198 376375 377062 377941 378139 378152 378289
## [3401] 378327 378436 379059 379913 380165 381393 381424 381477 381883 382122
## [3411] 382982 383044 383059 383209 383269 384206 384600 385034 386004 386764
## [3421] 387010 387722 388672 389739 389786 390008 390598 391290 391424 391872
## [3431] 392636 392814 392864 393286 394713 394716 394843 396019 397058 397628
## [3441] 397970 399216 399375 399612 401460 402544 402549 403834 404369 404900
## [3451] 405348 405475 405716

match(sub_cistrans_GoDMC_hvCpG_matched_control$controlCpG_name, cpg_names)

##    [1] 344658 182750 398917 293604 250517 389028 320321 361346 278675  79767
##   [11] 297360  20162 213019  28197  18957 355647 205134  51869 240454 163411
##   [21] 139238 255410 210252 153545 129311 253354 204799 341489 150143 171247
##   [31]  64826 308826 112906 119534 224191 316816 152619 405092 171396 231271
##   [41] 355301  10499 319793 361702 372847 355389 119559 217878  38334 215281
##   [51]  37565 332759 130027  44323 305624 322247 188489 206251  66157 225744
##   [61] 276414 323179 155969  78869 381722 354251 179738  95822  12968 242724
##   [71] 210861  59015 326099 266528 198514 236280  62355 164550 270759  28383
##   [81] 162353  25888  56635 190661  91778  37849 124291 302759 193166 131965
##   [91] 242476 238634 103454  52898 140345  86322 362635 403206 283927 196724
##  [101] 390594 272364 224530 142960 187232 102510 111968 273870 270559 266240
##  [111] 143657  54236 339136 158267  30176 405531 210207 383211 145098  82235
##  [121]  39541 360594  79602 139637 135957 196778 342448 223056 101168 374824
##  [131]  39245 205960 312703 290264 144169 394562 202013 357295   8565 133524
##  [141] 391589 301318  77075 266517 360856 297382 389037 256239 379427  12993
##  [151] 198871  69139 343335 123444 268052  81574   9334     21 300369 372471
##  [161]  71200 167982 328586  59677 349948 379155 378828 176975 184511 263079
##  [171] 156415   8903  86956 118211 196037  74284 209728 307900  37244 175669
##  [181]  89597  62676 101833 292163 374645 289260 167522  64375 156786 177712
##  [191] 377033 169584 382713 334161  71238 116008  33656  40916   9624 392015
##  [201]  44833 252486  39066 239606 156118 399966  78531 200276   1631  55311
##  [211] 127936  81509 278171  10767 318667  96679 199669 164779 351051 287063
##  [221] 142139 340134 231008 310644 255602 144496 278707 357860 296332 174829
##  [231] 268804 346802  83887 145670 299549 375166 290098 139185  26490 400222
##  [241]  60654 378768  59880 105148 288717 335927 129778 395765  21611     80
##  [251] 312473  17107  86810 395544 344231   4629 195331 400484  68568 141590
##  [261] 136102 357246 343365 215178 188908 340802 245503  78625     93  97795
##  [271]  38295  12448  48488  75216 279467  67965  90688 124808 109700  71164
##  [281] 343951 273866 205173  22168  29465 325200 359715  50891  65235 386050
##  [291]  79361 209992 151235 382820 163294 203927 398463 295297 366249 299635
##  [301] 112729 296451 367888 346820  24220 171299 315804 244847 234332  20867
##  [311] 204157 297626 131460  68844 347487  62664 204577   1462 372918  19269
##  [321]   3234 261513 344551 294342 296234 403949 106422 193189 227096 299277
##  [331] 282406 267134 171685 109784 374393 217122 110076 121867 376290 206737
##  [341] 151906 314896 397545   5356 136521 261555 252519 156412 302074 390179
##  [351] 145158 331949 153691 170539 154095  92232 216181 340418 221295  89733
##  [361]  59404 191500 166903  79923 369601 403097 282110 347634 107901 400062
##  [371] 291902 375713  32047 398526 352102 252154  75207 184244 217368  22962
##  [381] 401449 307993 175110 381067 281796 347929 323771   3528 382236 276781
##  [391] 261788 152560 178724  57069 143450 240705 386732 227770 333685 211499
##  [401]  56216 119230 234134 250337 222794 180202 218842 123955  31458 371432
##  [411] 165747 176348 323486 111488 299092  86806  10571  68641 212116 372094
##  [421]  24212  84646 395531 310179 175376 181850 169297 332767 288336 197938
##  [431] 270999 118448 263829 299045  50625 353494 324521  64882 353778 237207
##  [441]  36027 208009  31293 134491  86933 261463 271263 118060 390761 341884
##  [451] 213206 142779  48158 222732 183122 148214 139096  87476  33889 256256
##  [461] 268687 111397  52711 167125  51330 258579 374870 345819 354126 359361
##  [471] 258702 264615 263688 247446  21567 193124 386283 132846 224360  16264
##  [481] 133994 213665 212546  24422 300182 204292 390849 374201 320811  98165
##  [491]  10503 216199  94813 321920 102678 293116 185184 201055  75464 208530
##  [501]  97487  45237 191200 138041 342061 114395 184856 288778 111177 225221
##  [511] 157284 191991 293341 308112  51532 167974 331879 279547 305543 340486
##  [521] 145007 226829 260517 152990  48215  97326 315994 295782 342236  72977
##  [531]  44359 389640 311301 131650 307688 105600 167949 380540 303083 402317
##  [541] 284191 264001 346574  64635 102026 287094 146010 365194  74182 173485
##  [551]  31904 370379 148298 129921 111421 172288 350675 332934  38243 258742
##  [561]  40012  79854 369366 390506 369615 384436 239800  81267  54738  48972
##  [571] 113262 133019 296260 349951 239822 350517  96026 112770 101042  11214
##  [581] 342136 331195  16254 225540 181584 301397 233866 253554 222149 331109
##  [591] 227801 345393 306104 208766 317476  32151  48013 288792 362797 328843
##  [601] 174446 272928 212605  51244 269612 163574 404741 378564 375557 262230
##  [611] 298840 293364 375983 324340  72900 214141  97053 285143 273010 167036
##  [621]  79652 166504  69898 155820  53729 332222  82160  16049 171862  22525
##  [631] 214047 141837 192101 181637 137858 101007   3676 338481  55358 287464
##  [641]  26532 380893 379982  56103 205956 261290 301621 175036 157738 259517
##  [651] 225743 166265 399166 159152  23586 357778  45446 297095 122198  42861
##  [661] 200770  60247  28728 341243 254185  33437 170157 362275 272706  76548
##  [671] 157004 127558 169748 249834 247252 363309 283867 272940 215586 234226
##  [681] 233371 108713 398054 114625 289254 298707 238313 320499 405000 206585
##  [691] 183822 311448 126785 185129 379727  69933 381793  29226  33310 228500
##  [701] 166175 274465 193596 319718 391597 399330 329503 346132  60513 394956
##  [711] 141673 123807 262356 333146 214396  43892 228781 350175 252016  24300
##  [721]  65352   8073 206067 211633 239649 296348 253549 202870 321499 213511
##  [731] 193830 197275  79173 163910  52076 284766  89404  69137 142331 132342
##  [741] 217552 268988 107543 125985 165557 246188 172297 268697 158805 168142
##  [751] 289587  99763 246570 285152 382805 333281 174363  15686  31134 134480
##  [761] 391860 396884   2796 306466 224862  55657  72494 170035 372519 136138
##  [771]  37040 301647 131106 157237   2395 287259 212669 337739 383630 187828
##  [781] 168523 217296 101723 205788   2714  42650 185904 152381  85345 365345
##  [791] 286105 294469 192530 175026 168195 373075 343822 190414  22116 121547
##  [801] 143310 263549 146681 292839 262037  33690 266447 159873  10796 284491
##  [811] 274371 287840 130211 332536  40393 378882 384523 328209 333956  57899
##  [821] 360539 234398 388985  22291  30611 108113 337456 228311 125535 157744
##  [831] 156484  34964  38582  93513 104229 387259   8969 375416 343164   7995
##  [841]   9892 374690 307739 297630 376350 196512 281293 131508  95966 134537
##  [851] 253747 365964 236081 331417 333762 252022 310491 137646  20875 252008
##  [861] 227595 142675 127019  79259 116722  83415 372445  56888 101319  81891
##  [871] 249957 311762 218828 210267 170104 104685  54120 111000  41237  31277
##  [881] 266120  82367 278130 259411  86043 171951 263191   7140  85317  91532
##  [891] 173992  69179 365462 181398 197030 190700 240601 405854 110567  12395
##  [901]  92969 248476 403216 196363 164511 275912  24425  55157 138119 230759
##  [911] 108311 380422 215874 104490  31117 401651 303345 293771 117707 373765
##  [921] 280225 197087 288785 106393 356572 117427 358796 288725  31199 219329
##  [931] 102882 190878  11407  50049 197798  42821 285676 106109 338500 265655
##  [941] 218802 309686 246865   9285 255029 319338 185809  66387 300180 153994
##  [951]  60204 100755 252050 258077  88839 299080 262529 173450 268409 186149
##  [961] 285735 258925 163380 308135 182557 189102 325444 249888 174968  66674
##  [971] 137762 292571  57070 182741  60329 404246 371355 106272 169818 286516
##  [981]  12691 300927 385543 146582 362884 164360  21310 289911 230141  58122
##  [991] 140784 128051 251162  72726 351414 333437 106195 270956 200765   1936
## [1001] 253060 327422  76463  95301 142049  43831  36951 136839 398962  72816
## [1011] 375155 294123 244126 310565 242717 139310  67866  52529 114661  87832
## [1021]  89773 104477 126863  27969 388005  75875 147974 210952 136478  87687
## [1031]  82428 288848  20385 353511  68508  88633 358517 404572 351270  88803
## [1041]  92876 145149 402055 114924  37563  36349 244371 321463 228380 200102
## [1051] 265719  60749  98690 142911 369371  29475 312077 170881 218473  48899
## [1061] 198910 208068 155271  81423 136179 186619 176301  88484 358694  78450
## [1071] 403635 366301 274153 239331  57429 397566 318906 215728 131345 274782
## [1081] 179272 198016 259868  50641 295575 217435 187592 361844 119656 180642
## [1091] 219517 140626 278099  75417   7623 166330 339393 279007  54190 137979
## [1101] 115907    773 147382 403438 173773 367279 255802 238259 227946 351236
## [1111]  79529  91732 373167 311440 188103 232037 109922 159511  31375 403264
## [1121] 286258 251319 222173  33541  14464 190829 272828 244322 133743 138311
## [1131] 372843 309938 132526 286356 303332 276628  63419   5354 161450 308443
## [1141] 402791 163722  93040  76729 195018 178898 193822 182595 130936 179615
## [1151] 216903  91580 123440  92718  93634 240085   7484 117900 393872 339351
## [1161] 173083 365871 129341 214916  63306  48452 196031 321502 247220 208789
## [1171] 135516 150980 215802 298901 127768 127096  26524 168757 128680 363998
## [1181] 138473 143923 110483 367741 293543 386298 165417 185457  19242  55663
## [1191] 145581 357823 296900 308359  36051 278209 151820  17062 212505  79462
## [1201] 155302   4107 380975 103564 201014 198532 110012  15697 354180 125583
## [1211]  23861 282122  78118 307982 272267  83211 221708 170374 354337 216221
## [1221] 209564 114822 189924 244872 238506 196936 387537 255298  57526 305148
## [1231] 304510 392619 279354 253723 404817 174762  97059 344713 348601 355900
## [1241] 259866 254075 133540 133270  88160 138277 332830 256160  15843  35616
## [1251]  39057 287074 289946 151434 327091 118921  92397 232053 191770 311392
## [1261]  33064 330949   1503 261241  42470 200676 187620 285193   3371 162323
## [1271]  93881 382759 211478  34437 395282 113975 267228  30037 297524    527
## [1281]  57495 177782  15931 177412 325156 239896 222676 162076  95625 154902
## [1291] 335116 156180 303149  25880 207321 135543 366274 174076 272333  17323
## [1301] 141568 219685 240764 261066  92332 245357 345599 254659 309778 241730
## [1311] 363306 249832 165536 313466 245851 391029  45147 106502 354542 111497
## [1321]  55828  31062 324558 265659   6798 307110 184653   6598 106761  18797
## [1331] 173348   5883 366987 348680 159366 293980 351680  14406 357485 334205
## [1341] 203088 173441  77353  86832 183574 368429  14270 113373 333657   3938
## [1351] 275199 227206 241054 335353 193543 242483 246433 323704  22254  34316
## [1361] 158064 244857 345462 118821 288513  94014 111500 239367 118998  97460
## [1371]  20453 181855 250944 130378  18271  25116 290324  86035  43011 303514
## [1381] 303632 334001 287247 113888 180730 164042 224385 298041 358058 213916
## [1391] 360135 265182 231890 172097  11358 211719 304266 136858  60129 224016
## [1401]  90798 335629 228987 142409 201473 176799 369759 225294 318642 238234
## [1411]  24529 100853 377025 133924  74575  79153 197968   8202 267665  36562
## [1421] 224257  76535 157602 340253  24237  36576 330612  27448  54604 218656
## [1431] 378812  29044 361365 116215 347122 129929 348058 328073 104852 342250
## [1441] 222172 216119 126030 122160 218500  40983   9440 118750   2397 278180
## [1451] 171550 359604 299542 267140  26943 330776  25697 123096 320875 232260
## [1461]  55044 401364 380304 181274 220934  26353  27277 392625 111948 117145
## [1471]  77665 209043  11462 304769 343695 244474 355552 168412 262106  52018
## [1481]  87487 224435  67845 275216 261170 280073 296483 181291 269628 326777
## [1491]  90107 222806  21239 116284 326687 399466 291193 249236 123078 238732
## [1501]  99122 393648 158722 161125 280804 345707  63742 240493  83060 306258
## [1511]  11931   4920 265206 371541 151074 218307 310596  72965 348110  26986
## [1521] 174443   9495 206548 403326  82109 359426  53134 164658 356542 209262
## [1531] 317583 256315 207602 263755 240439 275674 346077 304522 234533 291544
## [1541] 337350  67983 120668 185272  65733 390148   5637 177451 382200 315840
## [1551] 345577 156747 203643  15299 349430 371990  16412 260586   9364 237261
## [1561] 251568 128719  36821  61942  23949 136918  27033 299789 338307 369208
## [1571] 244323 207390 189639 400053   1134 272882 311743 191415 152164 150438
## [1581]  28716 169378 219694 106700 340044 359962  87480 229858 125299 216158
## [1591]  34909 289530  99919 405813 402144 299314   4104 229508 302136  65513
## [1601]  92172 281982  48410 122207 380724 218340 160823 260822 166614 367119
## [1611] 240304 329907 145937 139557  37532 383667 307215 125087 128004 404790
## [1621] 180587 218787   2583 250115 106783  83154 169894 358608 341316 342823
## [1631] 244506 380024 277078  60385 106441 296833 146101 190466 277378 129830
## [1641] 216342 249356 125001 340088 252304 123177 225841 143874 371367 138090
## [1651] 302373 126656 159983   9561 121776 213126 325106  10387 101542 118105
## [1661] 186915   3494  97040 195578 300574 130551 246395 263475  18348 151776
## [1671] 284801  75237 390597  86056  56710 348631 276929 236267 139899  78127
## [1681]  32309  83871 392000 150621 304693 277004 277557   2846 230437 402635
## [1691] 289382 240094 116912 328189 185262 146956 321808   4209  19653 121690
## [1701] 113997 134252 175739  97895 100091  27321 103372  88729 207835 396125
## [1711] 119629 111403  81050  98384   1757 371252  37667 238510 100608  16569
## [1721] 300431 307271 274951 216531 328882  16562 172155 213638  13816 111123
## [1731] 129011 287772 265846 194709 187946 396980 136452 183277 361035  30762
## [1741] 210113  91116  56011 149226 170694 165465 140934 337618 199186  71805
## [1751]  35023 317917 229633 185605 265931 364647  88039  44454 130235 225656
## [1761] 164251   3697 389571 378938 243985 403111 291522  25840 336943 205758
## [1771]  58226  25042 400358 133981  27940 244121 388511 375725 342102  77347
## [1781] 133783 126025 193537 203006 148379 185961 164653  93424  38548 261519
## [1791]  24075 162640  13045 173695  38329 230192 117392 350401 237048 327999
## [1801]  17413 229017 392754  10751  13743  60146 396245 128850  33626 344790
## [1811] 121227  12569  48449 179924 291435  51702 232620 155684  33004 127553
## [1821] 358906 186615 252555 313797   5785 367694 117571 384414 197183 299797
## [1831]  65049 204149 327084 364110 346826  21410 178301 119550  65507  74120
## [1841]  78178 226531 245931 134856  16823 233113 357940 142226 394250 355240
## [1851] 166941  64747  33483 265464 110595 349774 214818 207442  65928 372353
## [1861] 116810 100423 146862 284638 176571 305696 382680 404685  55744 268364
## [1871] 403960  11375  98429 298286 105887 383565 314822 232684 272239 312602
## [1881] 334087 237842 240614 118063 246544 340237 107648 399509 387659  72586
## [1891] 187300 295980 113885 349947 236524 256540 138291 227912 212321 229855
## [1901] 203634   4748 155561 400523 353208 176637 171137 128173 143736 164993
## [1911] 146676 281301 329745 247569 286451 229056 217582  80993 231464  77670
## [1921] 341287 151635 199882 291624 259933 304295  82808 248184 150385  62789
## [1931]  70428 319222 394269 343576 358416  40087 225857 249964 332271 225876
## [1941]  76238 211263 154402 150609  42425 148695 147779 267333  97207 307710
## [1951] 255473  70342 118224 302738 300043 296086 267101  67440 118580 364336
## [1961] 295964 143447 399337 239209 215759 388952 359004 223487 396875 395476
## [1971] 119113 237675 176068 246992    501  55280 382528 127352 312097  32409
## [1981] 391931 355346 394526 258573 178642 194840 237628  39786 272714 271910
## [1991] 178643 258408 124929 186588 334442 114968  45387  19777 388508 101669
## [2001] 202887 324915 302934 340020 143311 119519 333588 194309 343533 174049
## [2011]  90393 101154 229136    401 301170 190432 230393 238569  59903 327593
## [2021] 166551 316035  42658  35528 112875 387853 158369  78829 118647  27494
## [2031] 186495 399100  39809 134977 221336 349923 298977 338252 376264 271048
## [2041]  58646 196133  93431 240109 225880 138585 222225 190885 239123 350299
## [2051]  79908 260847  39206 310248 241767 262301 394214  97693 106800 213777
## [2061]  85926  87680  79479  48012 296036 163484 354010 274809 293660  42373
## [2071] 131064  11971 208330 224389  47722 288808 159050 225230 340109 367107
## [2081] 230935 230074  53187 318537 265675 112701  89616 238090  40085   8407
## [2091] 176607  46783 187980 320480 130565  91040 326036 196289  11534 173988
## [2101] 285674  24665  84376 224506 224558 276595  96823  95716 333812 150586
## [2111] 185254 181313  56652 103338 194899  98618 199769   1924 383631 152881
## [2121]  80701 209650 163128   8789  60273 405044 338563 149254 335224 366923
## [2131] 244461 267218 118205  45294 199021 334815 373489  29762 243591 244235
## [2141]  39036 356621 259011 239178 248488 314365 150568  90892  29512 402520
## [2151] 330761 382149   3161 236718  85541  59525 323840 239164 363515 188508
## [2161] 396713 374579   5836 252156  57029 304900 117970 130585 345442 196336
## [2171] 116954 276682 156543 113395 277433 273875  34655 314652 371386  21839
## [2181] 208129 296901  81520  10846  21517 374838 100569 300493 121816 402898
## [2191] 129236 346824 376725 345753 222257 155000 380345 290964 388460  85512
## [2201] 308211 377893 396971 356959 397043  50901 304225  11932 274410 263269
## [2211] 375151  22174 316246  38459 339710 382223 304557 178178 377848 176491
## [2221] 152140 208931   7344 128478 231741  20715  36499 109855  55725  42360
## [2231] 190000  65356  37232 348521  95024 225574 382568 134503 127428 109534
## [2241] 205106 252552 279902 257197 349956 186865  72065 274572 104445 193335
## [2251] 405826 212375 387868 272220 293559 270461 279062 242458 178064 134385
## [2261] 350706 241358  83558 325181 167176 143093 297646 114660 297393 122145
## [2271] 244401 193882 220703 252691 367023 307708 150885 104768 315538 385632
## [2281] 181992 231124   3172  45301 373528 218068 165903   7769 184132 323854
## [2291]  17809  94265 207479 128288 371467    349 267483 154157 148618  43001
## [2301] 142573 170622 126062  70376 128886  99978 113744  29045 133715  73580
## [2311] 270186 265034 337495  51066 156149 381428 127761 144014 388739 359519
## [2321] 134327 286114 369842 357107 253418 203452 389094 213980 144211 228585
## [2331] 394221 161700  52946 396151 137871 266925 190553 321280 200209 124536
## [2341] 245488 101647 177699 306583  77898 151602  94865 114948 121389  57383
## [2351]  57847  65593 360086 226603 329213 366381 200885  32335 295562 224565
## [2361] 309501  95094  59537 391017  30633 156455 340388 366316 285929 137481
## [2371] 117221  17643   9445 104301 302562  13967  69561  98790 381789 102897
## [2381] 315488  23055 274408 365274 263369 348697  71492 391987 124015 172017
## [2391] 347894 304836 118907 233137   3073  10322 296829  47987 269328 257729
## [2401] 131771 119721  48043 331884  14222 302275 378544 366392 247060 232309
## [2411] 215003  26218 205633 204950 198372  27631  79381 348018 168227 322770
## [2421]  47403 213564  28692 323332  97883  85792 126880 312627 334898 200833
## [2431] 160301 305260 145318 142447  38113 244626 183040 374717 161451 156416
## [2441] 193010  42030 197519  49464 251598  48888 285570 322727  46095 301002
## [2451] 163249 132276 160355  48737 397063  16610 118615 140913 358205 114319
## [2461] 333606  73324 402725 399774 212661 392570 296556  95578 378863 219378
## [2471] 385278 289936 199258 197203 140935  47197 399858 336098  57315 376701
## [2481]  46289 281243 372514  99137  67921    847 168118 325148 368570  29124
## [2491]  72946 164814 396096  72281  71827  58611 331282 386265 271748  48633
## [2501] 272652  87193 269729  96684 287265  46510 368072 227254  96150 378769
## [2511]  59929 214321 357672  17958  22241 365184 196038 196979 278513 197570
## [2521]  53149 302773 241239 113636  18832 116438 150012 140766  59316  59043
## [2531] 357725 181769 341423 284043 370766 166363 371991 231607 183899  74961
## [2541]  50511 256604 323169 263396 237981  48762 357501  25300 173009 384550
## [2551] 107119 251793 127221 186955  87961  67628 248710  89592 373164 231578
## [2561] 188279 167780 358430 181635  81577 329963 372409  46260 382294 363531
## [2571]  86616 339946 101396 163874  30767 339560 161715  21026 393639 177889
## [2581] 133292 139659 138357 332516 360829 347360  47606 253918 214031 117109
## [2591]  69741  41632  91647 139253 202432  73110 103190 379605 387782 363378
## [2601] 271928 131945 178144 309989 146199 217566 273920  28116 172389  35373
## [2611] 122181 185992  71706 255033 114445  92951 161467 164892  50963 244178
## [2621]  23559 180086  49743 270129  56208 388136 392051 172470  19640  94416
## [2631] 184243   9401 367397 381940 116276 271784  15404   1268 387518 363312
## [2641] 255969  35520 313788   5767 211803 283292 331019  16453 226629 242219
## [2651]  32394 308721 210039 204408 260017 300663 169761 282825 173034  23936
## [2661] 392397 128070 351435  65045 397424 363104 391485 280195 373905  35900
## [2671]  10953 381398 297142 364078 125405  40227  40007  32788 289864  67612
## [2681] 312762 337288 118903 137616 230222 176323 171375 202484 337594 383507
## [2691] 169505  12027 301028 381859 104705  95087    335   3512 297125  74599
## [2701] 312109 162883  96048 403849 221366 400074 317022  66661  35838 126629
## [2711]   9331 387915 271643  60081 118865 287866 114234 209777 258973 293600
## [2721] 249598 235455 205416  59384  48357 362569 359286 403908 306396 201598
## [2731] 147129 335577  27129 120283 120006   7137 180454 171534  29252 295054
## [2741]  27174  96460 308067 359230 118085  87824  38668 244618 381244 191309
## [2751] 231476 102915   7542 290180 169989 253688 303072 208805 131575 168567
## [2761] 297754 176355 123782  31637 272669 348255  80189 208928  11679  62514
## [2771] 225136 304272 144816 305062  13491 254442 305306  60211 261810 163005
## [2781] 120260 391105 384546 211767  88593 152928 314309 262296 352255 313648
## [2791] 227423  38149 344631 143340 248223 122367 390111 125768 347702 149738
## [2801] 394007 185748 122043  64993 130765  64586 245397  69392 224509 308486
## [2811] 159486  78628  23360 369568 161130 167206  20597 116021  15264 158828
## [2821] 196209 262631  99181  66515 357581 295590 249047 137901 195738  65144
## [2831] 314560 237808  71033 309377 123750 197337 219171 312141  61370 307042
## [2841] 220858 353882  74000 259283 386207  33462 242348 340805 266909  73399
## [2851] 226855  97684 284778 118476 107495 226467  48041 315250 218449 183140
## [2861]  88047  65680 339520 262389  31876 290043 287434 380525 373309 192479
## [2871] 356814 301344  30316 405393 405334 151016 316112 260132 174343 118204
## [2881] 113211 265972 224776  81307  80859 245738  76356 399531 276867 373698
## [2891] 338107   5686 312423 390323  15810 322005  29146 304170  46841 382427
## [2901] 198048 315283 195417 123619 189721 321050 260194  34752 186046  47264
## [2911] 154562 115397  69522 196152  26946 310751   7694 314776  79251  32592
## [2921] 291876 154891  47017 381460   9516  75134  39690 304602 383920 281473
## [2931] 320628  62903 279850  63372   3598 299068  99926  28287 143692 250959
## [2941]  24997 346344 209039 179304 146082  85192 104153  40849  79977  73202
## [2951]  32291 202777 135412 349414 212293 171635 147217 296209 374015 255421
## [2961] 328304 206480  83750 317830 137079 239529 178535   7840 107033 131495
## [2971] 257830 375133 191395 296253 343354 369435 149808   4702  54633  70578
## [2981] 158755 157633 189917 136966 385394 100838 141962  37485 358581  39553
## [2991]  34493  84417 311569 130137 225761  74652  48010  34547  22296 362458
## [3001] 306973 235128 320293 288284 191382 374748 188433 137819 165554 385133
## [3011] 388039 215140 190337 236312 338153  44570 220517 179371 240571 308871
## [3021]  47845  61560 275383  75234 172397  85110 306173  42860 391879  15614
## [3031] 129494 347274 351409  29562 293822  64166 310875 114500 311347 372851
## [3041]   1737  80318 231631  77535 141116  70872  86917 372626 358672 219112
## [3051] 257117  43253 285016 311031 171244  11099  32869 236332 289963  46368
## [3061]  30537 221161 308449 399720  38128 337939 278083 214788  58655 116151
## [3071] 379687 309347 162878 264166 336295 173052 316174 134403  35806 275624
## [3081] 187977  84577  70019   5627   6779 376952 130694 308870 292851 336454
## [3091]  49801 184094 348161 190914  44739  95077  95936  53252 169126 252769
## [3101] 366822  71393  49636 273210 133754 236194 101107 336038 398999 387745
## [3111]  42640  97192 365611 137870 348884 203268 354451 278044 285043 354038
## [3121]  22426 103322 200653 269225 155163 256526 293784 245748 387671 204557
## [3131] 260345 240208 236613 138247 282133 107466 124753 275919 256599 121927
## [3141] 225866 131793 304530 223660  19648 326386 327323 258089 166858 221433
## [3151] 113418 354508  81417 335783  59841  23253  33978 176623 171661 309675
## [3161] 220598 393767 265945 405617 214364 387557  85642  90226 116927 215317
## [3171] 240182 390375  53680 125979 231997 204797  30063 282963 128071 301151
## [3181] 127020 216189 232753  39556  39548 260845 121782 146106 358641 262799
## [3191] 295963 153456 241260 345709 324421 360483 308204 311982 258651 259528
## [3201] 253817 316317  33206 208679  21503 262340 370091 317845  70222  70574
## [3211] 194909 322437 184032 315186 244713 146151 163238 351701 359152  73075
## [3221] 252763 182549  11779 204939  94762  17829 184536 281077 156148 404160
## [3231] 113187 384367   7223  67441 342942 129590  18799 325201 180307 312157
## [3241] 146262 135676 109407  41738 341575  58778 119632  87386 320200 154523
## [3251]   7620 313454  62965  40051 390462 377512 227201 329754 185615 224343
## [3261]  47835  39192 282391 220671  62981 104157 172730 329306 281056 108488
## [3271] 163453 127315 231097 126936  39194   9360 123339 180128 207994 197964
## [3281] 161562 298453   6527 215844  37105 342124 189890 293377 139682  11795
## [3291] 265088 358597  62170  15046  22179  25005 290034 124712 198987 180043
## [3301]  65597 327073   6561  94612  15454  75533 375706 342221 168826 284141
## [3311] 146889 380935 327699 196156 395572 208203 361136  76482 312755 274631
## [3321] 234558 310655  68260 110164 392908 370047 390608 369276 330686 258914
## [3331] 320595 322567  78418  77954 194307  12425 233377 259917 197141 282471
## [3341] 268838 126032 188906 138600  73190 129666  20670 187090 118911 372937
## [3351] 231557  92091 243789 183960 114910 187260 304866  86451 228431  43224
## [3361]  42929 329962 297158 337081 288172  68757 194889 350722 191065 193561
## [3371] 259462 155182 137746  52666 385383 110211 328582 190479  18477 352693
## [3381] 147120 144801 164481  11717 208652 219035 198925 124120  38465 259679
## [3391]  93778 300209 131085 304301  53034 222953  69733  15267 403589  71766
## [3401] 394365 404477 266053 262284 314787  43163 308559 118155 299500  54624
## [3411] 282191 212430 241723  77473 285828 220427 312767  33989 109152  59779
## [3421] 256448 349100  20280  91934  65458 141336 374904 184074  67151 179318
## [3431] 300462  73176 232672 216393  28300 114841  84963 281348 282600 187448
## [3441]  46444 329083 405550 381655 391266 128857 314420 305396 363119  48422
## [3451] 268897 248982 161879

pos2check <- c(match(sub_cistrans_GoDMC_hvCpG_matched_control$hvCpG_name, cpg_names),
               match(sub_cistrans_GoDMC_hvCpG_matched_control$controlCpG_name, cpg_names))

## Made with "Brent" and new method
load("/home/alice/2024_hvCpG/05_hvCpGalgorithm/resultsDir/Mariads/results_Maria_pos2check_0_8p0_0_65p1.RData")

Brent = data.frame(CpG = cpg_names[as.numeric(row.names(results_Maria_pos2check_0_8p0_0_65p1))],
                   alpha_Brent = as.vector(results_Maria_pos2check_0_8p0_0_65p1))

Brent %>% 
  mutate(ishvCpG = ifelse(CpG %in% Marias_hvCpGs$CpG, "hvCpG in Maria's study", "mQTL matching controls")) %>%
ggplot(aes(x = CpG, y = alpha_Brent, fill = ishvCpG)) +
    geom_point(pch=21, size = 3, alpha = .8)+ facet_grid(.~ishvCpG)

load("/home/alice/2024_hvCpG/05_hvCpGalgorithm/resultsDir/Mariads/results_Nelder_Mead_my_list_mat_Mariads_hvCpGandControls_0_8p0_0_4p1.RData")
NM = data.frame(CpG = row.names(results_Nelder_Mead_my_list_mat_Mariads_hvCpGandControls_0_8p0_0_4p1),
                alpha_NM_old = as.vector(results_Nelder_Mead_my_list_mat_Mariads_hvCpGandControls_0_8p0_0_4p1)) 

NM %>% 
  mutate(ishvCpG = ifelse(CpG %in% Marias_hvCpGs$CpG, "hvCpG in Maria's study", "mQTL matching controls")) %>%
ggplot(aes(x = CpG, y = alpha_NM_old, fill = ishvCpG)) +
    geom_point(pch=21, size = 3, alpha = .8)+ facet_grid(.~ishvCpG)

merge(Brent, NM) %>% 
  mutate(ishvCpG = ifelse(CpG %in% Marias_hvCpGs$CpG, "hvCpG in Maria's study", "mQTL matching controls")) %>% ggplot(aes(x=alpha_Brent, y=alpha_NM_old, fill = ishvCpG)) +
  geom_abline(slope = 1) +
  geom_point(pch=21, size = 3, alpha = .8) +
  theme_bw()

The new method is much better than the old one

ROC curves to define p0 and p1 matching Maria’s approach

\[ p0 = P(Z_{j,k}=0 \mid Z_j=0) \\ p1 = P(Z_{j,k}=1 \mid Z_j=1) \]

p1 = 0.65 in Maria’s approach (to be a hvCpG, it needs to be hypervariable in >65% of the datasets)

# Step 1: Read results and compute AUC + threshold
files <- list.files(
  "resultsDir/Mariads",
  pattern = "^results_Maria_pos2check_.*\\.RData$",
  full.names = TRUE
)

auc_df <- tibble(p0=double(), p1=double(), AUC=double(), threshold=double())
roc_list <- list()

for (f in files) {
  if (file.info(f)$size == 0) next
  nm <- load(f)
  res <- get(nm)
  
  pstr <- str_match(f, "_(\\d+(?:_\\d+)?)p0_(\\d+(?:_\\d+)?)p1")[, 2:3]
  p0 <- as.numeric(gsub("_", ".", pstr[1]))
  p1 <- as.numeric(gsub("_", ".", pstr[2]))
  label <- sprintf("p0=%.2f p1=%.2f", p0, p1)
  
  df <- as.data.frame(res) %>% rownames_to_column("CpGpos") %>%
    mutate(CpG=cpg_names[as.numeric(CpGpos)])
  colnames(df)[2] <- "alpha"
  
  truth <- ifelse(df$CpG %in% Marias_hvCpGs$CpG, 1, 0)
  roc_obj <- try(roc(truth, df$alpha, quiet=TRUE), silent=TRUE)
  if (inherits(roc_obj, "try-error")) next
  
  aucv <- as.numeric(auc(roc_obj))
  thr <- coords(roc_obj, "best", ret = "threshold", transpose = FALSE)[["threshold"]]
  auc_df <- auc_df %>% add_row(p0, p1, AUC=aucv, threshold=thr)
  
  pts <- coords(roc_obj, "all", ret=c("sensitivity","specificity"),
                transpose = FALSE) %>%
    as_tibble() %>%
    mutate(FPR = 1 - specificity, TPR = sensitivity, label=label)
  roc_list[[label]] <- pts
}

roc_all <- bind_rows(roc_list)

# Step 2: Plot ROC curves
roc_plot <- ggplot(roc_all, aes(FPR, TPR, group=label, color=label)) +
  geom_line(size=0.8, alpha=0.6) +
  geom_abline(slope=1, linetype="dashed", color="gray50") +
  labs(title="ROC Curves for Different (p0,p1)",
       x="False Positive Rate", y="True Positive Rate") +
  theme_minimal() + theme(legend.position="none")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Step 3: AUC heatmap with thresholds
best <- auc_df %>% slice_max(AUC, n=1)

heatmap_plot <- ggplot(auc_df, aes(p0, p1, fill=AUC)) +
  geom_tile(color="white") +
  geom_text(aes(label=round(threshold,2)),
            size=3, color="black", check_overlap = TRUE) +
  # Outline best tile
  geom_tile(data = best, aes(p0, p1), fill = NA, color = "red", linewidth = 1.5) +
  scale_fill_viridis_c(option = "plasma") +
  labs(title="AUC & Optimal Threshold per (p0,p1)",
       x="p0", y="p1") +
  theme_minimal()
# Step 4: Show both plots
print(roc_plot)

print(heatmap_plot)

The best threshold ploted in tiles is the alpha probability that gives the best sensitivity-specificity trade-off. It’s derived using a statistical rule (Youden’s J: maximizes (Sensitivity + Specificity – 1), giving the best overall balance between true positives and true negatives). It classifies a CpG as hypervariable in Maria’s approach if alpha > threshold

Conclusions (NB TO UPDATE WITH NEW ALGO):

• Best AUC = for p0 (true negative) is quite low in Maria’s data (0.45)
• Threshold for this best AUC = : Classify a CpG as hvCpG in Maria’s data if alpha ≥
• Previous discrepancies between methods are likely due to scaling!

Subset Maria’s data following Atlas data structure (less N/dataset) to check power of detection in Atlas data:

3 tries for (p0, p1):

• c(0.45, 0.6): best AUC, closest to Maria’s data; low specificity and mid sensitivity * c(0.95, 0.65): our original idea, high specificity but mid sensitivity * c(0.9, 0.9): high specificity and sensitivity

makeMyPlots <-function(p0p1){
  df1 = get(paste0("results_Nelder_Mead_my_list_mat_Mariads_test3000CpGsvec_", p0p1))
  df2 = get(paste0("results_Nelder_Mead_my_list_mat_Mariads_mimicAtlas_test3000CpGsvec_", p0p1))
  
  p_str = str_match(p0p1, "(\\d+_\\d+)p0_(\\d+_\\d+)p1")[, 2:3]
  p0 = gsub("_", ".", p_str[1])
  p1 = gsub("_", ".", p_str[2])
  label = paste0("p0=", p0, ", p1=", p1)
  
  plot1 = data.frame(df1, df2) %>%
    tibble::rownames_to_column("CpG") %>%
    mutate(ishvCpG = ifelse(CpG %in% MariasCpGs$CpG, "1500 hvCpG in Maria's study", "1500 mQTL matching controls")) %>%
    dplyr::rename(alpha_normalDS = alpha_hat, alpha_mimicAtlas = alpha_hat.1) %>% 
    ggplot(aes(x = alpha_normalDS, y = alpha_mimicAtlas)) +
    geom_point(aes(fill=ishvCpG), pch =21, size =3, alpha = .3) + 
    geom_abline(slope = 1, linetype = "dashed") + theme_bw() +
    ggtitle(label = label)
  
  # === 1) Prepare the ORIGINAL ===
  df_original = df1 %>%
    data.frame() %>%
    tibble::rownames_to_column("CpG") %>%
    mutate(ishvCpG = ifelse(
      CpG %in% MariasCpGs$CpG,
      "1500 hvCpG in Maria's study",
      "1500 mQTL matching controls"
    ),
    source = "original") %>%
    melt()
  
  # === 2) Prepare the MIMIC ===
  df_mimic = df2 %>%
    data.frame() %>%
    tibble::rownames_to_column("CpG") %>%
    mutate(ishvCpG = ifelse(
      CpG %in% MariasCpGs$CpG,
      "1500 hvCpG in Maria's study",
      "1500 mQTL matching controls"
    ),
    source = "mimicAtlas") %>%
    melt()
  
  # === 3) Combine ===
  df_combined = bind_rows(df_original, df_mimic)
  
  # === 4) Plot ===
  plot2 = ggplot(df_combined, aes(x = value, y = ishvCpG, fill = ..x..)) +
    geom_density_ridges_gradient(scale = 3, rel_min_height = .01) +
    facet_wrap(~source, ncol = 1) +
    scale_fill_viridis(name = "alpha", option = "C") +
    theme_bw() +
    theme(axis.title.y = element_blank()) +
    labs(title = "Probability alpha of being hypervariable in my algorithm",
         subtitle = "Original data of Maria vs mimic Atlas structure")
  
  print(plot1)
  print(plot2)
}

1. low specificity and mid sensitivity (p0=0.45 and p1=0.6)

best AUC, my algorithm with the closest proximity to Maria’s approach)

# load("resultsDir/Mariads/withAlpha0.5atStart/results_Nelder_Mead_my_list_mat_Mariads_mimicAtlas_test3000CpGsvec_0_45p0_0_6p1.RData")
# makeMyPlots(p0p1 = "0_45p0_0_6p1")

2. high specificity and mid sensitivity (p0=0.95 and p1=0.65)

high true negative rate, i.e. high specificity)

# load("resultsDir/results_Nelder_Mead_my_list_mat_Mariads_mimicAtlas_test3000CpGsvec_0_95p0_0_65p1.RData")
# makeMyPlots(p0p1 = "0_95p0_0_65p1")

3. high specificity and sensitivity (p0=0.9 and p1=0.9)

# load("resultsDir/results_Nelder_Mead_my_list_mat_Mariads_mimicAtlas_test3000CpGsvec_0_9p0_0_9p1.RData")
# makeMyPlots(p0p1 = "0_9p0_0_9p1")

Notes:

lambdas are defined based on a 5% threshold, even if alpha is not. How to not hardcode them? MCMC on lambda? Link with alpha?

The analysis is complete.